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ABSTRACT 

We develop a numerical scheme and code for estimating the energy and mo- 
mentum transfer via neutrino pair annihilation [v + v ^ e~ + e"*"), bearing in 
mind the application to the coUapsar models of gamma-ray bursts (GRBs). To 
calculate the neutrino flux illuminated from the accretion disk, we perform a ray- 
tracing calculation in the framework of special relativity. The numerical accuracy 
of the developed code is certificated by several tests, in which we show compar- 
isons with the corresponding analytical solutions. Using hydrodynamical data in 
our coUapsar simulation, we estimate the annihilation rates in a post-processing 
manner. We show that the neutrino energy deposition and momentum transfers 
are strongest near the inner edge of the accretion disk. The beaming effects of 
special relativity are found to change the annihilation rates by several factors in 
the polar funnel region. After the accretion disk settles into a stationary state 
(tj^jically later than ~ 9 s from the onset of gravitational collapse), we find that 
the neutrino- heating timescale in the vicinity of the polar funnel (< 80 km) can 
become shorter than the hydrodynamical timescale, indicating that the neutrino- 
heated outflows can be launched there. We point out that the momentum transfer 
can play as important role as the energy deposition for the efficient acceleration 
of neutrino-driven outflows. Our results suggest that the neutrino pair annihila- 
tion has a potential importance equal to the conventional magnetohydrodynamic 
mechanism for igniting the GRB fireballs. 
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Introduction 



Accumulating observational evidences have been reported so far, identifying a mas- 
sive stellar collapse as the origins of long-duration gamma-ray bursts (GRBs), such as the 
preponderance of short-lived massive star formation f or their host galaxies, as well as the 



identification of SN Ib/c light curves i n the afterglows (jPaczynskil Il998l : iGalama et al.lll998 



Hjorth et al. 


2003; 


Stanek et al. 


2003) 



the accretion of debris falling into the central black hole (BH) (jPiro et al.l 119981 ) . Pushed by 



(Wooslev 


1993; 


Paczvnski 


1998; 



In the collapsar model, the central cores with significant angular momentum collapse 
into a BH. Neutrinos emitted from the accret i on di s k heat the matter of th e polar funnel 
region to launch the GRB outflows. iPaczynskil fjl990l ): iMeszaros fc ReesI (119921 ) pioneerlingly 
proposed that the energy deposition proceeds predominantly via neutrino and antineutrino 
annihilation into electron and positron (e.g., u + u — )■ e~ + e"*", hereafter "neutrino pair 
annihilation"). The relativistic outflows are expe cted to ultiin ately form a fireball, which is 
good for explaining the observed afterglow (e.g., iPiranI ( 1l999l )). In addition, it is suggested 
that the strong magnetic fields in the cores of order of 10^^ G play also an active role both for 
driving the magnet o-driven jets and for e xtracting a significant a r nount of energy from the 
centra l engine (e.g.. IWheeler et al.l ( 120001 ); [Thompson et al.l (l2004j ) ; lUzdensky fc MacFadyen 
(120071 ) and see references therein). 



However, it is still controversial whether the generation of the relativistic outflows pro- 
ceeds predominantly via magnetohydrodynamic (MHD) or neut rino-heating p rocesses. S o 



far, m u ch attention has been paid to the MHD processe s fe.g.. jProga 



roga et al. 



(2003); Mizuno et al.l (120061) ; iFuiimoto et al.l (120061 l2007h ; iNagataki et al.l (l2007f ); iNagataki 
(120091 ); iHarikae et al.l ( 120091 )). A general outcome of these extensive MHD simulations is 
that the magneto-driven shock waves, produced mainly by the field-wrapping of the strong 
precollapse magnetic fields or by the field amplification due to magnetorotational instability 
(MRI), can blow up massive stars along the rotational axis. Although those primary jet-like 
explosions are at most mildly relativistic due to too much baryons , they may be related 



to the X-ray flashes, which is a low energy analogue of the GRBs (ISoderberg et al.l 12006 



Ghisellini et al.l 120071 ). and those jets, which make a low density funnel along the rotational 
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axis, could provide a good birthplace for the subsequent relativistic outflows ( iHarikae et al. 
2009h . 



In contrast to such blossoms in MHD studies, there are only a few studies pursuing the 
possibility of generating jets by the energy deposition via neutrino pair annihilation. This 
is mainly because the neutrino emission from the accretion disk generally becomes highly 
aspherical, thus demanding one in principle, to solve the mul tidimensional neutrino t r ansfer 
problem. For the first time in the numerical simulations, iMacFadyen fc Woosleyl (119991 ) 
pointed out the importance of the energy deposition via neutrino pair annihilation in their 
simulations, however the energy deposition rates to the polar funnel region were adjusted 
by hand to produce jets. Without those artificial energy injections, neutrino-driven outflows 
have been yet to be realized so far in the numerical simulations, to our best knowledge. It 
should be noted in the MHD-driven coUapsars that too much precollapse magnetic fields 
with rapid rot ation can lead t o the MHD-driven explosion, wh i ch can prohibit t he for mation 



of a BH (e.g., IPessart et all (120081 ): iBucciantini et all ( 120091 ): iBurrows et all (120071 )). thus 



not good for the collapsar scenario. Therefore it is still important to address the possibility 
of the neutrino-driven mechanism in collapsars. 

Thus far, several methods for the implementation of the neutrino pair annihilation, 
have been reported. One straightforward method is to solve the full radiation field by the 
Monte-Carlo method (lTubbslll978l : iJanka fc Hillebrandtlll989l ). Still at present, this seems 
too expensive to combine with the hydrodynamic simulations. By estimating t he local fluxes 



and sp ectra of the neutrino emission via the so-called neutrino leakage scheme, iRuffert et al. 



(119971 ) proposed to estimate the heating rate by summing up the local contributions of the 
neutrino and antineutrino radiation incident from all directions. However this scheme re- 
quires the double angular integration for neutrino and antineutrino for all the grid points in 
three dimension, wh ich is still highly t i me-co nsuming and prevents its application to hydrody- 
namical simulation. iRuffert fc Jankal ( 119981 ) have upgraded this scheme to be more feasible, 
by defining the position of the neutrinosphere. I nstead of the three- dimensional stellar vol- 
ume as neutrino and antineutrino sources like in IRuffert et al.l (Il997l ) , the energy deposition 
rates can be estimated by summing only over the two-dimensional neutrinospheres, which 
reduces the computational cost significantly. 



Along this prescription, iNagataki et al.l (120071 ) have estimated the neutrino heating 
rates, and included them to the hydrodynamical simulation. For reducing the computa- 
tional time, they adde d one more assumption of the optically thinness of the accretion disk 
to the prescription by IRuffert fc Jankal (119981 ). Even with this potential overestimation of 
the heating rates, th e neutrino-driven outflo ws were not observed in their simulations. As for 
the methodology of IRuffert fc Jankal (119981 ). the neutrino fluxes are estimated by summing 
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up the local rates only along the z— direction (: perpendicular to the equatorial plane), which 
simplifies the directional dependence of the incoming neutrinos. The directional dependence 
can be more appropriately handled by doing the so-called ray-tracing analysis. Moreover, 
even when the spacial structure of the accretion disk becomes highly inhomogeneous, which 
is the case for our collapsar simulations, the ray-tracing analysis is good for omitting the 
surface of the accretion disk, from which the neutrino rays cannot travel, such as from the 
regions opposite to the central black hole. 



More recently iDessart et al.l (120091) have developed a new scheme to estimate the energy 
deposition rates v ia neutrino pair annihilation using the mutliangle neutrino-transport solver 
f lOtt et al.ll2008al ). They discussed the possibility of the neutrino-driven outfiow production 
in the postmerger phase of binary neutron-star coalescence. Without the general relativistic 
effects, this is the state-of-the-art method to estimate the neutrino pair annihilation. As 
for the general relativistic effects, iBirkl et al.l (120071 ) have performed a series of ray-tracing 
calculations in a Kerr spacetime. They found that the total energy deposition rates mea- 
sured by an observer at infinity can increase by a factor of two by the general relativistic 



effect, which is consistent with the result of previous studies (e.g.. lAsano fc Fukuyamal ( 12001 



Joining in the extensive efforts mentioned above, here we present a numerical scheme 
and code for calculating the neutrino heating via pair neutrino annihilation, which is de- 
signed to be incorporable to th e hydrodynarn i cal si m ulations of coUap s ars. R elying on the 
neutrino leakage scheme as in iRuffert et al.l ( 119971 ): iRuffert fc Jankal ( 119981 ). we estimate 
the incident neutrino flux by doing a ray-tracing calculation from the neutrino sphere, on 
which the neutrino distribution function is assumed to take a Fermi- Dirac distribu t ion. T he 
ray-tracing calculation is done in the Minkowskian spacetime (e.g.. iKotake et al.l (j2009aj)). 



neglecting the general relativistic (GR) effects (e.g.. iBirkl et al.l ( 120071 )) for simplicity. Thus 
our scheme is limited to capture the special relativistic effects such as beaming effects. One 
of the main characteristics of our method is that the neutrino pair annihilation can be esti- 
mated only with quantities deflned on the numerical grids of the hydro dynamical simulation. 
This can be done by carrying out a Lorentz transformation of radiation variables in the co- 
moving frame back to the Eulerian frame. By doing so, the computational cost of the double 
angular integration that is required in estimating the annihilation rates, can be reduced sig- 
niflcantly. We derive the formulation required for this, and describe it elaborately. We check 
the numerical accuracy of the developed code by showing the comparisons with analytical 
solutions, some of w hich we newly give in this paper. Based on the results of our long-term 
collapsar simulation (IHarikae et al.ll2009l ). we run our new code to estimate the heating rate 
in a post-process ing manner. Although the presen t ed sch eme is not the state-of-the-art in 
comparison with iBirkl et al.l (120071 ): iDessart et al.l ( 120091 ). we elaborately study the possi- 



Fig. 1. — An example of the ray-tracing (white hne) of neutrinos for estimating the neutrino 
pair annihilation towards a given point outside the accretion disk (indicated by red). The 
central sphere represents the black hole (BH). We trace each neutrino trajectory backwards 
till it hits to the surface of the neutrino sphere, which closely coincides with the surface of 
the accretion disk (footpoints of the rays on the torus). Note again that the geodesic of 
neutrinos is given by the straight line, since the general relativistic effects are not included 
in this study. 

bility of the neutrino-driven mechanism in collapsars and will point out that the neutrino 
pair annihilation has a potential importance equal to the conventional MHD mechanism for 
igniting the GRB fireballs. 

This paper is organized as follows. We summarize the special relativistic formulation of 
the neutrino energy deposition and momentum transfer in Section [2l Section 3 is devoted to 
the main results. The numerical tests and its comparison with the analytical solutions are 
shown in section 3.1. In section I3l2l we estimate the annihilation rates in a post-processing 
manner using hydrodynamical data in our collapsar simulation. We summarize our results 
and discuss their implications in Section HI 
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Fig. 2. — Schematic picture of the two coordinate system bridging between one for the 
hydrodynamic quantities (O) and another for the radiation quantities (O'). Note that the 
coordinate O' is used to describe the target position, where the neutrino pair annihilation 
occurs. denotes the neutrino/antineutrino source. rs,u{'rs,v,ds,v,4>s,u) is the position of 
Su measured in O, while ri^{r^, 9^, 4>u) is the position of 5*,^ measured in O'. r(r, 9, cj)) is the 
position of O' measured in O. 

2. Neutrino pair annihilation in special relativity 

Before going into details, we first show an illustrative figure, which presents our calcu- 
lation concept for the neutrino pair annihilation in collapsars (Figure [1]). For estimating the 
neutrino pair annihilation at a given target, we trace each neutrino trajectory (white lines) 
backwards till it hits to the surface of the accretion disk (colored by red). It should be noted 
that the surface closely coincides with the surface of the neutrino sphere (footpoints of the 
rays on the tori), from which neutrinos and antineutrinos emerge freely out. As indicated in 
this figure, some fraction of neutrinos comes out from the the inner edge of the disk, whose 
rotational speed is close to the speed of light. This suggests that the special relativistic 
effects become important towards the precise determination of the heating rates. 

For later convenience, we define the two frames namely of O and O' in the spherical 
coordinate systems (Figure [2]). Here O is the coordinate system that we use for describing 
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the variables obtained in the hydrodjTiamical calculations. On the other hand, O' is the one 
for describing the target points where the pair annihilation occurs. Note that the coordinates 
between O and O' system, differ only in the position of their basis. 5*,^ in the figure represents 
the position of neutrino/antineutrino source. is usually described by the O' system. In 
the following, variables measured in the O' system, are described with the subscript u or z/, 
indicating that they are related to the neutrino sources. With these definitions, we move on 
to describe the special relativistic formulation from the next section. 



2.1. Formulation of energy deposition and momentum transfer rates 



The energy dep o sition rate via neutrino pair annihilation (e.g., iGoodman et al.l ( 119871 ): 



Asano fc Fukuyamal (120001 )) is written as 



dtdV 



fu{,Pv,r)f^{py,r)a\vy - Vp\{e^ + ep)d p^d Pp, (1' 



where is the number density of neutrinos in phase space, p^, e^, v^, is the momentum, 
energy, and velocity of neutrino, respectively. Those definitions are the same for antineutrino 
by changing the notation u to u. a is the cross section given by 

a = 2c'KGl{p,-p,), (2) 

where the dimensionless parameter K is written as 

l + 4sin2ew + 8sin^^w 
K (z/e, i/e) = , (3) 

DTT 

_ _ 1 - 4 sin^ 6'w + 8 sin^ 6'w , 
K{iy^,, Vf,) = K{ur, Vr) = • (4) 

DTT 

Here, the Fermi constant Gp = 5.29 x 10~^^cm^ MeV~^, and the Weinberg angle is sin^ 9^/ = 
0.23. Note in our setting of the coordinate system (e.g.. Figure 2) that r corresponds to the 
origin of the O' (targets) system measured in the O system. Variables with subscripts of u 
(z/) are quantities upon the neutrinosphere, which is labeled by Si,(p) in Figure 2. They are 
also measured in the O' system. 

Aiming at a straightforward incorporation of the ray-tracing calculation into the hy- 
drodynamical simulations, we calculate the energy deposition(momentum transfer) rates in 
the laboratory (lab) frame. On the other hand, the quantities related to radiation, such as 
the distribution function and emissivity /absorptivity are locally defined in the rest frame 
of fluid, in which the radiation isotropy is maintained. Designating the subscript to the 
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variables in the rest frame , the Lorentz transformations between the two frames are given in 



Mihalas fc MihalasI fll984f ) as, 



dtdV 


= dtodVo 


de 


= —deo 






dVt 


= 4^^o 


e 


= 7 (1 + /3yUo) eo 




= eo/b{l-(3fi)], 



(5) 
(6) 

(7) 



where /3 = |V|/c, 7 = (1 — and fi = V ■ n/lVl with V and n, being the velocity of 

fluid and the unit vector describing a given ray of neutrino, within a solid angle of dQ. 

Substituting from Equations ([5]) to ([H]) in Equation (1) gives the heating rate in the lab 
frame as, 

^^^j^ = 2cKGl I d9,dcP,d9,d(P, 
dtdV J 

X [1 — sin 9u sin Op cos {(f^ — fu) — cos O^, cos sin 9^, sin 9p. (9) 
Here reflects the special relativistic correction to the neutrino energy as 

= l/[7,(l-/i,/3,)]. (10) 

Here it should be noted again that the variables with the subscript and u denote the ones 
defined in the rest frame of the neutrino source and on the neutrino sphere respectively, such 
that 



(3u = (11) 
c 

is the velocity of the fluid on the neutrino sphere with its Lorentz factor of 
and the directional cosine of 

^ \V I ' ^ ^ 
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where n,y is the direction of the rays of neutrino from the neutrino sphere, 



\Pv\ 



(sin 6'i, COS 01/, sin 6'i, sin </);/, cos 6^1, 



(14) 



The following two quantities are the energy-weighted integration of the neutrino distribution 
function on the neutrino sphere (namely fvfliXufi^Pvo)) ^^s. 



4,oUo{ru,o,P^fl)de 



(15) 
(16) 



As in lRuffert et al.l (119971 ) ; iRuffert fc Jankal ( 1l998l ) ; iNagataki et al.l (120071 ) , the neutrino num- 
ber flux is simply assumed to be conserved as fu{,f)Pv) = fyi.f'u,o,Pyo) along its trajectory 
to the target. The neutrino distribution function on the neutrino sphere is assumed to take 
a Fermi-Dirac shape with a vanishing chemical potential as. 



U{ru,o,Pu,o) 



1 



driQ 



(hc)^ deodilodtodVo 
1 1 



(17) 



where T^^ is taken to be the same with the matter temperature on the neutrino sphere of 
T{ru). Now the energy integrals in Equations (llSp and (ITB]) can be expressed by the Fermi 
integrals J-^ as 



My) 



X 



exp(x — y) + 1 



dx. 



{he 
{hcf 



-^3(0). 



(19) 
(20) 



The surface of the neutrino spheres is defined per each lateral angular grid in the hydrody- 
namic calculation [9) as 

/ -dr = 2/3, (21) 

where R^^ {9) is the radius of the neutrinosphere for the neutrino species of i, namely for z/e, i>e, 
and i'ji.T {pfi.,v) with the corr e spond ir ig mean free path calculated by a multi-flavor l eakag e 
scheme JEpstein fc Pethickl Jl98lh : iRosswog fc Liebendorfeil J2003h : iKotake et aP J2003h : 
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Takiwaki et al.l (|2009[ )) with special relativistic corrections ( iHarikae et al.ll2009[ ). As for the 



neutrino opacities, electron capture on proton and free nuclei, positron capture on neutron, 
neutrino scattering with nucleq n and nuclei, photo-pair, plasma processes are included (e.g.. 



Rosswog fc Liebendorferl (120031 )). 



Following the same procedure above, the momentum transfer rate in special relativity 
can be also derived as, 



drn}. 



dtdV 



2KGl 



de,. 



^udOpC 



X [e^(r, n,)^t{r, n,)E,,o{r, fi.)iVp,o(r, fip)n, 
+e^(r, n,)C{r, fi.)i?p,o(r, l^p)iV.,o(r, fi.)n,] 

X [1 — sin 9u sin 9p cos {ip^, — (fp) — cos 6i, cos 9p]^ sin 9i, sin 9p. (22) 

Finally the remaining procedure is to change the coordinates of the angle integrations in 
Equations ([9]) and ( l22l) from the O' to O coordinate system. Albeit sacrificing the formulae 
to be messy (e.g., appendix A), this procedure, by which the neutrino source and the target 
can be connected one by one, merits for reducing the numerical costs significantly. The 
coordinate transformation with respect to the solid angle is written as, 

~l" Jfitpif^i^^ 'Pv^ f^s,uj 'Ps.u) \rs,u=ri^pi^ ,^dflg ud(j)g u 



+ Jiprifl'U,' 



(23) 



where r^,^(rs,,., 



^u) is the position of 5",^ measured in the O system. Using the expres- 
sions of Jacobian (Jr^, J^^, J^^, see appendix A for details), we can calculate the annihilation 
rates based only on the quantities defined on the hydrodynamical grid-points. 



2.2. Implementation of ray-tracing calculation 

We have shown how to calculate the annihilation rates when neutrinos from the neutrino 
sphere can reach directly to the target region. However in reality, neutrinos can be absorbed 
in some situations before hitting to the target, for example, by encountering with the optically 
thick region. To correctly count the real contribution, we set the following three criteria by 
utilizing the ray-tracing calculation. 

The first task is to define an outward neutrino emission from the neutrino sphere. This 
can be satisfied hy by ■ dy > 0, where the former is the normalized neutrino momentum 
from the neutrino sphere dy{= Py/\Py\ = (r — rg^y)/\r — Vg^yl), and the latter is the normal 
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direction to the differential surface of the neutrino sphere. As mentioned, we define the 
surface of the neutrino sphere to be dependent only on the lateral direction (e.g., Equation 
f ET]) ). Thus the analytic form of b^, can be simply given using the valuables such as Og^i, and 
(f)s,u that denote the surface position of the neutrino sphere. Furthermore t he criterion of 



by ■ dy > would provide more realistic regulation of the neutrino flux than in iRuffert et al. 



(119971 ): iRuffert fc Jankal (119981 ) whose criterion is determined only by the density gradient 



at the neutrino sphere. 

The second criterion is the positivity of (2/3 — f) along each trajectory from the surface 
of the neutrino sphere to the target region. Here f denotes the (locally-defined) opacity 
for each neutrino species (e.g., the integrand of the left-hand-side of Equation (|2T]) ). This 
is necessary to omit the neutrino flux which comes from the optically thick region, such as 
inside the neutrino sphere. For not counting radiation absorbing to the black hole (BH), 
the third criterion is the positivity of r — Vg, where Vg = 2GMbh/c^ with G, Mbh being the 
gravitational constant and the mass of the BH, respectively. For simplicity, we will set Vg 
to be the Schwartzshild radius when we apply the ray-tracing calculation to the coUapsar 
simulations in the later section. However in reality, the BH is rotating so that the accurate 
estimation of the geodesies is necessary, which can be an another possible extension of this 
study. 



3. Result 

3.1. Numerical Tests 

Before applying the newly developed code to collapsars, we shall check the accuracy of 
our code. For the purpose, we derive some analytical solutions for test problems and present 
comparison with the numerical results in this section. 

Amongst a list of tests, the flrst requirement we set is to check whether the newly 
derived Jacobians (Equation (23)) and their numerical implementations are correct. This 
will be checked by reproducing radiation flelds shedding from a light-bulb in spherical and 
non-spherical geometry (section 13.1.11 and 13.1.21) . The second one is to check the accuracy 
of our ray-tracing calculation. Note that this can be checked in each test, because the 
ray-tracing is done in whichever tests. In addition, we present an additional test, which 
is specially designed to test the validity for the application to the coUapsar simulations in 
section 13.1.31 

In the following tests, our numerical domain has the spherical coordinates with equally 
spaced 100 radial meshes covering from to 500 km in radius. To check the numerical con- 
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vergence, we vary the numerical resolution for Q and </) direction, respectively. For simplicity, 
the axial and equatorial symmetry is assumed for the hydrodynamic quantities (such as for 
density, electron fraction and so on) as a background. The density and temperature on the 
neutrino sphere is set to be p = 10^^ g cm~^ and T = 5 MeV, respectively. 

To measure the deviation from the corresponding analytical solutions, we estimate the 
error function of i?, which is defined by the Li norm as. 



E 


~ r ' 


(24) 






-^l,err 


— ^ 1 -^num 3; ana | ; 


(25) 


-^l,ana 


= S^^ana) 


(26) 



where Xnum represents the numerical values such as the energy deposition rate (gjp) and 
the absolute value of the momentum transfer rate (|m^p|), while Xana represents the corre- 
sponding analytical solutions. S indicates the summation, which is taken over the whole 
computational domain. 

As for the momentum transfer rate, the error function of E cannot always be a useful 
measure because the analytical solution becomes zero in several situations (see sections 13.1.21 
and I3.1.3|) . In such a case, we define another error function of El as. 

El = (27) 

-'"max 

Mn,^ = S|m+|, (28) 
M^ax = Sm+,,, (29) 

where ?7imax is the analytical solution, which is defined from the energy deposition rate as 
'^max = ?a!iia/C; thus meaning the maximum momentum transfer. With this £'2, we will 
evaluate how precisely the cancellation can be maintained. 
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Fig. 3. — Energy deposition (left panel) and momentum transfer rates (right panel) from 
a spherical neutrino sphere (central white circle). In each panel, left-half and right -half 
represents the numerical and analytical solutions, respectively. The numerical resolution is 
set to be {ne,n^) = (60,64), which are indicated in the top right edge in each panel as 
"60X64". With Figure 4, our code is shown to reproduce the sphericity of the radiation 
fields well. 
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Fig. 4. — Deviation from the analytical solutions for the energy (left) and momentum transfer 
rates (right) as a function of and ng in the case of the spherical neutrino sphere. See 
Equation fl2^ for the definition of error which is calculated by the Li norm. 



3.1.1. Radiation from Spherical Neutrino Sphere 

As a simplest test, we first present how accurately we can reproduce the radiation field 
from a spherical source, which emits radiation isotropically from its surface. In this test, 
we set a spherical inner boundary of 50 km in radius, inside which material is taken to be 
completely optically thick for neutrinos, while outside is taken to be completely optically 
thin. Albeit very crude, such a situation is akin to the neutrino radiation from the neutrino 
sphere in the postbounce phase of core-collapse supernovae. In addition to this static source, 
we consider another case in which the spherical neutrino sphere rotates relativistically. By 
this test, we hope to see and test the special relativistic effects, which is one of the main 
focus in this paper. 



Non Relativistic Case 



In the case of the static neutrino sphere, t he analytic formu l ae of t he heating and momentum 
transfer rate have bee r i expl i citly given in iGoodman et al.l (119871 ) ; ICooperstein et al.l ( 119871 ) ; 
Salmonson &: Wilson! Jl999l ): lAsano fc Fukuvamal J200oh 



as 



dqUr) AcKGl {kTf 



dtdV 



{hey 



■J-4(0)J-3(0)F(r) 



(30) 



and 



(im+p(r) _ AKGl {kTf 
dtdV ~ 3 {hcf 



■J-4(0)J-3(0)G(r)n, 



(31) 
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where geometrical factor F{r) and G{r) is given by 

27r2 



F(r) 



{1 - Xy{5 + AX + X' 



and 



with 



TT 



Gir) = —(1 - X)\l + X)(8 + 9X + SX^) 
o 



X= Wl 



where Ri, is the radius of the neutrino sphere. 



Ru 



(32) 
(33) 
(34) 



Figure |3] shows the energy (left panel) and momentum (right panel) rates numerically 
calculated with our code (left-half) and the corresponding analytical solutions from Equation 
( 130|) (right-half). Note that the label such as "20X16" at the top right edge indicates the 
resolution of our ray-tracing calculation, here ng = 60, = 64, the equally spaced mesh 
numbers in the lateral (6) and azimuthal direction (0), respectively. Comparing the two 
panel, we find no visible differences in each panel. In Figure HJ the deviation from the 
analytical solution for the different numerical resolutions (e.g.. Equation (IMl) ) is shown. 
As shown, the relative error decreases with np^ns as low as .1 %. Compared t o oth er 



studies treating the radiation problems (e.g., IStone et all fll992[ ): [Turner fc Stond (120011 )) 



the obtained accuracy here seems high enough. These results assure the validity both of the 
newly derived expression of J^^ (from Equation f lAip to flA16l) ) and its implementation. 



Relativistic Case 

For the neutrino sources with relativistic motion, it is not generally possible to write down 
the special relativistic factor of C,u (Equation (10)) in an analytic form. This is because it 
depends locally on the angle between the velocity field and the flight direction of neutrinos. 
To see special relativistic effects, we examine a special case, in which the analytic solution 
can be found. The annihilation rates along the polar axis, emitted from the relativistically 
rotating neutrino source, can be analytically given as, 

d,Mr) _ icKGl (kTf ^^(0)^3(0)^,,), ,35) 



dtdV 3 79(/ic)6 

dmUr) AKGl {kTf 

-Hdv- = ^?(M^-^^^°^-^^^°)^(^)"- 

where the geometrical factors F{r) and G{r) are just same as Equations fl32|) and fl33|) . 
respectively. The sharp suppression with a factor of 1/7^ is the outcome of the relativistic 
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Fig. 5. — Energy deposition rate (left) and radial momentum transfer rate (right) from the 
relativistically rotating neutrino sphere (central white circle). Comparing with Figure [3], the 
energy and momentum rates are shown to be significantly reduced along the rotational axis by 
the beaming effect, (t is noted that the numerical resolution taken here is {rig, n^) = (60, 64). 



beaming, simply because the rotational velocity is perpendicular to the polar direction. We 
will see it soon in the following. 

Figure shows the energy deposition (left) and momentum transfer rates (right) from 
the relativistically rotating neutrino sphere. Here we assume 7 = 2 with vanishing Vr and 
Vg. Note in this figure that the numerical results are only shown, owing to the inaccessibility 
of the analytical solution as mentioned above. It is shown that the annihilation rates both 
for energy and momentum is greatly reduced near on the rotational axis as expected from 
Equations fl5B]) and fl36p . while they are enhanced around the equatorial plane. Figure 
E] shows the comparison of the numerical solution (solid line) with the analytical solution 
(dashed line) along the rotational axis for energy (left) and momentum (right) transfer rates. 
In support of the validity of our special relativistic implementation, no visible differences are 
seen between the analytic and numerical solutions. Figure [7] shows the deviation of E (e.g.. 
Equation (1241) ) for the energy (left) and momentum (right) from the analytic solution. From 
this test, it is shown that the lateral grid points {uq) should be at least more than 20 to 
maintain the percent levels of agreement with the analytic solution. 
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Fig. 6. — Comparison between the numerical (solid line) and analytic (dashed line) solution of 
the energy (left) and momentum transfer rates (right) for the case of the relativistic rotating 
neutrino sphere. It is noted that the numerical resolution taken here is {n0,n^) = (60,64). 
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Fig. 7. — Same as Figure H] but for the case of the relativistic rotating neutrino sphere. 



3.1.2. Spherical Cavity 

In this test, we set a spherical cavity of 300 km in radius, outside which is the uniform 
opaque neutrino source. Although this is not a realistic astrophysical situation, it provides 
one of the most simplest null test of the momentum transfer because the neutrino flux is 
isotropic everywhere inside the cavity. For the heating rate, this test merits because the 
geometrical factor (Equation fl30l) ) is analytically given by 

647r2 

F{r) = -—. (37) 



Comparing the numerical (left-half) and analytical solution (right-half) in each panel 
of Figure [8], the discrepancy is shown to become smaller for the better numerical resolution 
(from left to right). From Figure [9] showing the El and E2 errors for the energy (left) and 
momentum transfer rate (right), it can be seen that the errors mainly depend on n^. By 
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Fig. 8. — Same as Figure [3] but for the energy deposition rates in the case of the spherical 
neutrino cavity. 




Fig. 9. — Same as Figure H] but for the test of the spherical neutrino cavity. The error of E2 
(right) is calculated with the Li norm, which is defined in Equation fl29|) . 

setting n<^ = 32, we can get the agreement with an accuracy of 1 % for E and 0.1 % for E2, 
which would satisfy the requirement of the vanishing momentum transfer sufficiently. 

3.1.3. Non-Spherical Cavity 

We now move on to the non-spherical cavity test, by which we can check if the radiation 
contribution from the face element of drdcj) is accurately estimated. For simplicity, we just 
modify the shape of the spherical cavity in the previous section. We set the position of the 
neutrino sphere to change with the polar angle 9. For 6 < ir/A and 6 > Sn/A, the neutrino 
sphere is situated at 300 km, while for n/A < 9 < 37r/4, at 100 km. Region inside the 
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Fig. 10. — Same as Figure [3] but for the energy deposition rates in the case of the non- 
spherical neutrino cavity. 



neutrino sphere is taken to be completely thin, while completely thick outside. Even with 
the geometrical difference, it should be noted that the items to be checked are the same for 
the case of the spherical cavity, because the net flux is isotropic and does not change. 

From Figure [TOl it can be seen that some inhomogenities seen closely at the edge of 
the boundary (left-half: numerical) become smaller for the higher resolutions (to the right 
panel) . From Figure [111 the decreasing rate of the relative errors for E and E2 is shown to 
be smaller than the ones from the previous tests due to the relatively complicated geometry. 
By taking the resolution of {rig, n^) > (32, 32), our code can reproduce the analytic solutions 
with percent levels of accuracy here. 

As evident from the numerical convergence with increasing the numerical resolution 
(e.g., see Figures 4, 7, 9, and 11), the previous tests have supported the validity of the 
newly derived expression of the Jacobians (from eg. (lAip to (1A16P ). its implementation, 
and the ray-tracing calculation. From those test calculations, it is suggested that the grid 
numbers of (n^, "n.^) > (30, 30) are required for obtaining the percent levels of agreement with 
the analytical solution. Considering the computational cost, we choose to employ the grid 
numbers of (n^, uq, n^) = (300, 40, 32) in the actual implementation for coUapsar simulations 
in the following section. 
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Fig. 11. — Same as Figure |9] but for the case of the non-spherical neutrino cavity. 
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3.2. Application to Collapsar 

Having demonstrated the accuracy of our code in the previous sections, we now proceed 
to show an application to collapsars. In this section, we estimate the annihilation rates 
in a post-processing manner based on our simulation results of the long-term evolution of 
collapsar and then discuss the possible impacts on the collapsar dynamics. 



3.2.1. Recipes for Collapsar Simulation 



For obtaining the hydrodynamic profiles for the ray-tracing calculation, we perform 
the special relativistic simulations of collapsar ( iHarikae et al.ll2009l ). Without repeating the 
numerical procedures again, we will shortly summarize major modifications added in this 
work, namely about the initial model and the position of the inner-boundary. 

Since our focus is not on the MHD-driven collapsars, we construct initial models with no 
magnetic fields. Taking the precol l apse density /electron-fraction/entropy profiles of progen- 
itor of 350C in lWoosley &: Hegerl (120061 ) . we embed the angular velocity r2(r, 6) analytically 

as 

r]o^o' + «^iso(M(X))X4 



Q{r,e) 



(38) 



where X = rsin6', fiiso(M(X)) = jx^^^M {X)) / X"^ . M{X) is the mass coordinate at X, and 
jiso and r^iso is the specific angular momentum and angular velocity in the last stable orbit. 
Model parameters are a, fig, and Xq. This distribution provides co-rotation with = fig 
ai X < Xq connected to constant specific angular momentum j = aj\so at X > Xq. Xq 
would correspond to the size of co-rotating regi on in pre-coUapse phase , such as Fe core. 
To be closer to the original rotational profile in IWoosley fc Hegerl (120061 ). we set a = 0.8, 
Qq = 1.2 rad s~^, and Xq = Rpc, where Rpe is the size of Fe core (~ 3000 km for 350C). 

Another major modification is the position of the absorbing inner boundary of the 
computational domain. For reducing the computational cost, t he radius of the bo undary, 
which mimics the surface of black hole, was set to be 50 km (IHarikae et al.ll2009[ ). This 
manipulation may lead to the underestimation of the neutrino luminosity by excising the 
radiation from more inner edge of the accretion disk. This is obviously disadvantageous for 
producing the neutrino-driven outflow. In this paper, we take more compact inner-boundary 
as max(10 km, 2rg), extending to the outer boundary of 30000 km. 



- 22 - 



-2e+07 -le+O? 1e+07 2e+07 -IbtWl -le+07 lc+07 2e+07 




X[cm] X[cm] 



Fig. 12. — Typical hydrodynamic configuration when the accretion disk is in a stationary 
state (here, 9.1 s after the onset of gravitational collapse). In the left panel, the logarithmic 
density (in g cm~^, left-half) and temperature (in K, right-half) are shown. The right panel 
shows the electron fraction (left-half) and the Lorentz factor (right-half). The velocity, 
indicating by white arrows, is normalized by the scale in the bottom right edge of the box 
(here, 3 x 10^" cm s~^). The white solid line denotes the area where the density is equal to 
10^^ g cm""^, representing the surface of the accretion disk. The central black circle (11.5 
km in radius (~ 2rg)) represents the inner boundary of our computations. 



3.2.2. When can the pair neutrino annihilation work ? 



For energizing the jets via pair neutrino annihilation, the higher neutrino luminosity 
from the accretion disk and the lower density along the polar region, are favorable. Our 
coUapsar model starts with a rapid mass accretion into the central objects. Within ~ 2.0 s 
after the onset of gravitational collapse, this model accretes more than 2~ ?)Mq in the center. 
This indicates the black hole formation in the center because the maximum mass of the 
neutro n star, which the Shen e quation of state employed here can sustain, is in the same mass 
range (IKiuchi fc Kotakell2008l ). Simultaneous with the rapid infall, the neutrino luminosities 
also show a drastic increase, and then shift to a gradual increase reflecting the mass accretion 
to the newly formed accretion disk. At ~ 8 s from the onset of gravitational collapse, the 
total neutrino luminosities gets as high as 10^^ erg s~^, which consists of 2.2 x 10^^ erg s~^ 
for z/g; 9.4 X 10^^ erg s~^ for i>e, 5.1 x 10^^ erg s~^ for and Vt-. Afterwards, the total 
luminosity keeps almost constant with time, reflecting that the accretion disk comes to a 
(quasi)stationary state. 
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Fig. 13. — Same as Figure 12 but for the energy and momentum transfer rates with (left) 
and without special relativistic corrections (right). In each panel, logarithmic contour of the 
energy deposition rate (left-half) and the momentum transfer rate multiplied by the speed 
of light (right-half) both in unit of (erg s~^ cm~^) are drawn. 



Figure [12] shows the distributions of density, temperature, electron fraction, and Lorentz 
factor of our model at 9.1 s after the onset of gravitational collapse. From the left panel, the 
structure of the accretion disk with a polar funnel along the rotation axis is shown. The disk 
temperature reaches typically ~ 10 MeV, which provides the high neutrino luminosity of 
-|^g53 gj,g g-i_ pYom the right panel, the Lorentz factor at the inner edge of the disk becomes 
7 ^ 1.3, indicating the importance of the special relativistic effects. Taking these hydrody- 
namical data, we estimate the neutrino pair annihilation by the ray-tracing calculation in a 
post-processing manner. 



3.2.3. Energy deposition and momentum transfer rates in collapsar 

The left- and right- hand sides of the left panel of Figure [13] shows the energy deposition 
and the momentum transfer rates (its absolute value) . The energy and momentum transport 
from the accretion disk (regions colored by black) to the polar funnel regions (regions colored 
by red) can be apparently seen. Interestingly, the transfer rates become highest near the 
inner edge of the accretion disk. Inside ~ 80 km in radius from the center, the typical energy 
deposition rates along the polar axis are 10^° erg s~^ cm~^, where the momentum transfer 
rates are the order of 10^^ erg cm~^. Thus the momentum transfer rates multiplied by the 
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speed of light c becomes comparable to the energy deposition rate there. 



For our model, the density along the polar funnel regions is as low as 10^~^ g cm""^ 
(see Figure [TWleft)). Assuming that the momentum transfer there is maintained for 10 ms 
with conserving its momentum, we can estimate that the material along the funnel can be 
accelerated up to the Lorentz factor of 10^~^. This suggests that the momentum transfer 
plays as important role as the energy deposition for the efficient acceleration of the neutrino- 
driven outflows. The right panel of Figure [T3] is the annihilation rates without the special 
relativistic correction [^j^ = 1 in Equation (fTOl) ). Comparing the two panels in Figure [T3l we 
flnd that the special relativistic effect decreases the transfer rates by several factors inside 
the polar funnel. As already discussed in section [3. l.lj this is the outcome of the relativistic 
beaming by the rapidly rotating accretion disk. On the other hand, the reduction in the 
total deposition rates, given by the volume integral of the local rates in Figure [I3l can stay 
relatively smaller, namely of 6.22 and 7.62/(10^*^ erg s~^) with and without the corr ection, 
respectively. These values are comparable to what expected by lNagataki et al.l (120031 ) . since 
the mass accretion rate in this epoch is about O.IMq s~^. The resulting conversion efficiency 
of the heating, which is defined by the ratio of the total deposition rates to the total neutrino 
luminosity, is 0.514% and 0.630 % with and witho ut the correction. T h ose numbers are 
comparable t o the ones in t h e pre v ious studies (e.g . . iRuffert et al.l ( 119971 ): iRuffert &: Janka 



( Il998l . ll999r ): iNaravan et al.l ( 120011 ): ISetiawan et al.l ( l2004l . l2006l )). an d this negat i ve eff ect of 
special relativity on the energy deposition rate is also mentioned by iBirkl et al.l ( 120071 ). 



3.2.4- Comparison of times cales 



Based on the annihilation rates in the last section, we compare various timescales in 
this section, such as the neutrino-heating timescale and the dynamical timescale. Then we 
anticipate if the neutrino-heating outflows could be produced or not in the polar funnel 
regions. 

To trigger the neutrino-heating explosion, the neutrino-heating timescale should be 
smaller than the advection timescale, which is characterised by the free-fall timescale in 
the polar funnel regions. This condition is akin to the condi tion of the su ccessful neutrino- 
driven explo sion in the case of c ore-collapse supernovae (e.g.. iBethel (Il990l ) and see collective 
references in I Janka et al.l ( l2007l )). The heating timescale is the timescale for a fluid to absorb 
the comparable energy to its gravitational binding energy to be gravitationally unbound, 
which is deflned as Thcat = P^/q^, where p is the average density at a certain radius and 
we take p{r) = 3M{r)/A7Tr^, $ is the local gravitational potential. The dynamical timescale 

depicts the ratio of the dynamical t^j^ to the 



is deflned as Tdyn gure 
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heating timescales Theat, showing that the ratio becomes greater than unity inside 80 km in 
the vicinity of the rotational axis. This indicates the possible formation of the neutrino- 
driven outflows, if coupled to the coUapsar's hydrodynamics. 

Figure [15] shows various timescales for the material along the rotational axis in Figure 
UM Trci = pc^ /q^ characterizes the timescale for matter to become relativistic by the neu- 
trino heating. As mentioned, the first criterion of Tdyn > Theat, is satisfied for the regions 
inside ~ 80 km (compare green and blue lines). More interestingly, inside ~ 50 km, Trei 
gets slightly shorter than Tdyn (compare red and blue lines), indicating that matter could 
attain relativistic motions there. To see these outcomes needs the implementation of the 
ray-tracing calculation in the hydrodynamic calculation, which we pose as a next task of 
this study. It is noted that with our choice of the grid numbers ((n^, rig, n^) = (300, 40, 32)), 
it generally takes several minutes (in the CPU time) for each ray-tracing calculation, us- 
ing 256 processors in the Cray XT4 system at the National Astronomical Observatory in 
Japan. To follow the dynamics typically later than 10 s after the onset of core-collapse for 
satisfying the Tdyn > Thcat condition, 10000000 hydro-steps are generally needed because the 
hydrodynamical timescale is an order of 10"^ s in our simulation. Therefore it is still compu- 
tationally expensive to turn on the ray-tracing calculation throughout the whole simulation. 
In the actual implementation, we plan to switch on the ray-tracing calculation by monitor- 
ing intermittently Tdyn/Theat during the simulations in the computational domain. After the 
ray-tracing calculation is turned on, seeing Tdyn/Thcat ^ 1, the neutrino-driven outflows will 
be launched within 0.1 s (Harikae et al. in preparation). Thus the computational time using 
the above supercomputer can be reduced to several months, which are expected to be quite 
feasible. 



4. Summary and discussion 

Bearing in mind the application to the coUapsar models of gamma-ray bursts, we have 
developed a numerical scheme and code for calculating the energy and momentum transfer 
via neutrino pair annihilation in the framework of special relativity. To estimate the transfer 
rates, we perform a ray-tracing calculation of the neutrino flux from the accretion disk 
of collapsars. We have tested the numerical accuracy of the developed code by several 
numerical tests, in which the comparisons with the corresponding analytical solutions were 
shown. Using hydrodynamical data obtained in our coUapsar simulation, we have run our 
code to estimate the annihilation rates in a post-processing manner. We have obtained the 
following results. 

For the progenitor chosen in this study, the accretion disk settles into a stationary state 
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typically ~ 9 s after the onset of gravitational collapse. The accretion disk later on can 
be luminous source of neutrinos as high as ~ 10^^ erg s~^. At this epoch, the polar funnel 
along the rotation axis is formed, which is heated intensively from the disk via neutrino pair 
annihilation. The beaming effects of special relativity are found to decrease the transfer 
rates by several factors in the polar funnel regions. Inside ~ 80 km in the vicinity of the 
polar funnel, we point out that the neutrino-heating timescale can become shorter than 
the hydrodynamical timescale. This indicates the possible formation of the neutrino-heated 
outflows there. Our results also suggest that the momentum transfer via neutrino pair 
annihilation plays as important role as the energy deposition for the efficient acceleration of 
the neutrino-driven outflows. 

Finally we shall make comparison with previous work and also mention some imperfec- 
tions of this study. As mentioned, general relativistic effects ignored in this study have 
been considered as an impo r tant ingredient of the neu t rino pair annihilation in c ollap- 



sars (IJaroszvnskil Il993l . Il996l : ISalmonson fc WilsonI Il999l : lAsano fc Fukuyamal 1200 ll . 12000 



Birkl et al.ll2007l ). It is noted tha t we have ca l culate d the energy and momentum deposition 
rates for the thin disk models in iBirkl et al.l (120071 ) to assess the validity of our code. Our 
code can reproduce their results well both for the Newtonian and special relativistic case (see 
Figure [T6l for the 2D conflgurations) . In fact, the total heating rates are = 2.0 x 10^^ and 
1.5 X 10^®(erg s~^) f or the Newtonian (model DN) and the special relativistic case (model 



DN4L), w hich are in Birkl et al.l (120071 ). 1.5 x 10"^*^ and 1.3 x 10'^^(erg s~M, respect i velv ( see 
Table 1 in Isirkl et al.l J2007h l When we boldly extrapolate results in Isirkl et al.l J2007h to 
ours, the local heating rates can be enhanced significantly (by several factors) in the vicin- 
ity of the central BH due to the GR effects. To update the present scheme to the one in 
general relativity requires not only the ray-tracing in the curved space, but also the general 
relativistic formulation of radiative transport, which we are to investigate as extension of 
this study. At the same time, we should improve the ray-tracing calculation to take into 
account the charged-current absorption occurring along a given ray, which is ignored in this 
study. Without the neutrino absorption, the heating rates estimated in this paper should 
give a upper bound. We also plan to include the effects of magnetic fields on this study. 
By changing the strength of magnetic fields systematically, we hope to clearly understand 
how the outfiow production in collapsars, depending on the precollapse rotation, change 
from the neutrino-driven mechanism to the MHD-driven one. The magnetic effects on the 
annihilation rates a nd its impacts on the dynamics are also remained to be studied (e.g., 
Zhang fc Da] toO^ ). 



Furthermore, the neut rino oscillation by Mikheyev-Smirno v -Wolf enstein (MSW) effect 
(see collective references in Kotake et al. ( 2006 ): Kawagoe et al. ( 20091 )) could be important, 
albeit in much later phase than we considered in this paper. It should be remembered that 
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the typical density in the polar funnel was p ~ 10^ g cm~^ in the present model. When the 
density there drops as low as p < 10'^ g cm~'^ later, the neutrino oscillation could operate 
for neutrinos traveling from the accretion disk to the polar funnel. If this is the case, the 
incoming neutrino spectra to the polar funnel regions and the pair annihilation rates there 
could be af f ected significantly. Together with the effects of neutrino self-interaction (e.g.. 



Duan et al.1 (120061 ) ). there seem a lot of issues to be clarified about the neutri no oscillation 



i n the c o llapsar env i ronmen ts. As in the case of core-collapse supernovae (e.g.. iKotake et al. 



(l2009bl ): lOtt et al.l (l2008bl ) and see references therein), studies of gravitational- wave emis- 
sions from coUap s ars rn i ght provide us a new w indow to probe into the central engine (e.g.. 



Hiramatsu et al.l (120051 ): ISuwa fc Murasd (120091 )). As a sequel of this work, we are planning 



to implement the ray-tracing calculation in the hydrodynamic simulation and clarify those 
aspects. We just stand at a starting line to study those interesting topics, which will be 
presented elsewhere one by one soon. 
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Appendix A. Calculations of Jacobian 

In this section, we give an exact expression of the Jacobians ( J^^, J^</,, J^r) namely. 



drs 



^^ls 



dfls.iy d<j>s,v 



dii3,v dtps, 



(Al) 

(A2) 
(A3) 



by which we can estimate the heating/momentum-transfer rates only by the quantities de- 
fined on the numerical grids of the hydrodynamical simulations. 



From a geometrical calculation, albeit a little bit messy, the following partial derivatives 
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can be found straightforward as, 
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where the following valuables are dependent only on the position of the neutrino sphere 
{9s,v, 4>s,u) and the direction of a specified direction {9, 0) as, 
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Finally, it should be noted that above equations can be very simple along the polar axis as, 
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where Z is the distance from the equatorial plane. 
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Fig. 14. — Same as Figure [13] but for the logarithmic contour of the dynamical timescale 
Tdyn divided by the heating timescale Thcat- This is at the epoch of 9.1 s after the onset of 
gravitational collapse, when the neutrino luminosity is as high as 10^^ erg s~^. 
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Fig. 15. — Plots of various timescales for material along the rotational axis in Figure 14 (see 
text for the definition of the timescales) . 
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Fig. 16. — Anniliilation rates for the thin disk models in iBirkl et aLl (120071 ) (for model DN 
(left) with the Newtonian ray-tracing and for model DN4L with the special relativistic one 
(right). The value of the energy deposition rate (Q^, e.g., Equation (9)) is color-coded, while 
the spatial vector is visualized by showing the spatial velocity vector v = Q^^ jQ^ (Q^^ is 



equal to Equation (22)), which is normalized by the speed of l ight (c = 1) 



jeing represented 



by the arrow (top right in each panel). A good agreement with lBirkl et al.l (120071 ) is obtained 
by comparing to their Figure 2 (top right). 



